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EXECUTIVE  SUMMARY 


This  report  describes  the  basic  physics  and  numerical  techniques  used  in  implementing  a 
range*dependent,  tropospheric  microwave  propagation  model.  This  propagation  model  al¬ 
lows  for  both  finite  surface  cond^ictivity  and  variable  (i.e.,  highly  irregular)  surface  teiTain. 
The  model  is  based  upon  the  electromagnetic  parabolic  wave  equation  (PE)  and  employs 
a  novel  implementation  of  the  split-step  Fourier  PE  (SSFPE)  algorithm  to  efficiently  com¬ 
pute  the  electromagnetic  radiation  fields.  Tliis  physics  report  reviews  the  fundamentals 
of  the  PE  method  and  provides  a  detailed  error  analysis  of  the  SSFPE  algorithm.  The 
generalization  of  the  SSFPE  algorithm  to  handle  nonfiat  boundaries  is  also  discussed. 

The  propagation  model  is  implemented  as  the  VTRPE  (variable  terrain  radio  parabolic 
equation)  computer  code  and  is  used  in  the  prediction  of  radar  system  performance.  The 
VTRPE  code  has  the  following  characteristics: 

(1)  full- wave  propagation  physics  (i.e.,  field  amplitude  and  phase  are  computed); 

(2)  direct  solution  of  electromagnetic  fields; 

(3)  exact  treatment  of  refraction  and  diffraction  phenomena; 

(4)  exact  treatment  of  multipath  phenomena; 

(5)  range-dependent  atmospheric  refractivity  inputs,  N(z,r); 

(6)  infinite  or  finite  conductivity  surface  boimdary  conditions; 

(7)  linear  transmitter  field  polarization  (vertical  or  horizontal); 

(8)  variable  surface  terrain  elevation  and  surface  dielectric  properties; 

(9)  frequency  dependent  atmospheric  attenuation; 

(10)  frequency  range:  ss  0.1  -+  30  GHz; 

(11)  generalized  transmitter  radiation  patterns; 

(12)  arbitrary  transmitter/receiver  geometry; 

(13)  automatic  selection  of  the  SSFPE  range  step-size  and  FFT  transform  size;  and 

(14)  automatic  monitoring  of  SSFPE  solution  global  error. 

The  VTRPE  model  properly  accounts  for  the  dominant  mechanisms  governing  tro¬ 
pospheric  radio  wave  propagation,  including  the  effects  of  anomalous  propagation  arising 
from  spatial  changes  in  atmospheric  refractivity  and  variable  terrain  features. 
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§1.0  INTRODUCTION 


This  report  describes  the  physics  and  numerical  techniques  used  to  model  microwave 
propagation  in  range-dependent  environments  by  means  of  the  parabolic  wave  equation 

(PE). 

The  method  of  parabolic  wave  equations  was  first  proposed  in  1944  by  Leontovich^  a.s  a 
means  of  solving  elliptic  partial  differential  wave  equations.  He  used  the  technique  to  solve 
the  problem  of  electromagnetic  wave  propagation  above  a  plane  earth.  In  1946  Leontovich 
and  Fock^  applied  the  PE  method  to  the  problem  of  transhorizon  radio  wave  propagation 
above  a  spherical  earth,  thereby  making  a  breakthrough  in  electromagnetic  wave  propaga¬ 
tion  modeling.  Approximately  30  years  passed  before  a  practical  algorithm  for  solving  the 
Leontovich-Fock  parabolic  wave  equation  was  developed.  In  1973,  Hardin  and  Tapper!^ 
developed  the  split-step  Fourier  parabolic  equaton  (SSFPE)  algorithm  and  applied  it  to  the 
problem  of  modeling  ionospheric  radar  propagation.'^  The  split-step  Fourier  PE  algorithm 
exploited  advances  in  computer  hardware  and  the  development  of  the  fast  Fourier  trans¬ 
form  algorithm  to  yield  an  efficient  numerical  solution  to  the  Leontovich-Fock  parabolic 
wave  equation.  In  1977  Tappert^  introduced  SSFPE  methods  to  the  underwater  acoustics 
community,  where  it  rapidly  became  a  valuable  tool  for  predicting  range-dependent  under¬ 
water  sound  propagation.  In  1983  the  split-step  Fourier  PE  method  was  reapplied  to  radar 
propagation  by  Ko,  Sari,  and  Skura®  to  study  cuiomalous  tropospheric  microwave  prop¬ 
agation.  Subsequently,  Dockery  and  Konstanzer'  applied  the  SSFPE  method  to  analyze 
phased  eirray  radar  performance.  Since  then  several  others,  including  Ryan®  and  Craig, ^ 
have  developed  electromagnetic  PE  models. 

The  scope  of  this  report  is  restricted  to  a  description  of  the  electromagnetic  PE  model 
physics  and  the  split-step  Fourier  PE  algorithm.  A  detailed  description  of  the  numerical 
and  computationed  issues  associated  with  implementing  specific  algorithms  is  left  for  a 
later  document.  By  intent,  only  the  salient  and  most  important  ideas  are  described; 
no  attempt  is  made  to  cover  classical  electromagnetism  or  the  mathematical  theory  of 
partial  differential  equations.  The  reader  is  directed  to  standard  texts  for  the  necessary 
background  material. 

This  report  focuses  on  the  propagation  of  radio  waves  in  an  inhomogeneous  atmosphere 
with  the  treinsmitter  and  receiver  located  near  the  surface  of  the  earth,  in  the  troj)osphere. 
The  emphasis  will  be  on  short  wavelengths,  the  so-called  microwaves,  which  are  charac¬ 
teristic  of  radars  and  short-range  microwave  communication  links.  The  wavelength  range 
of  interest  is  from  about  3  meters  down  to  about  1  cm,  which  corresponds  to  a  frequency 
range  of  approximately  100  MHz  to  30  GHz.  This  encompasses,  to  a  large  extent,  the  bulk 
of  the  search  and  acquisition  radars  in  current  use.  It  will  be  assumed  that  the  wavelength 
is  sufficiently  small  that  ionospheric  reflections  are  absent  but  yet  long  enough  that  large 
numbers  of  atomic  or  molecular  resonances  in  the  gaseous  components  of  the  atmosphere 
do  not  occur  in  a  small  wav'clcngth  inter’ll.  In  other  w'ords,  it-  will  be  assunred  that  the 
frequency  is  well  above  the  ionospheric  plasma  frequency  but  low  enough  that  dispersion 
effects  are  not  important. 

A  further  assumption  is  that  the  atmosphere  can  be  modeled  as  a  simple  material 
medium,  viz.,  a  linear,  isotropic,  nonionized  medium.  This  implies  a  neglect  of  the  effects 
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of  the  eaxth’s  magnetic  field  on  wave  propagation,  and  a  restriction  of  propagation  paths 
to  the  lower  portion  of  the  atmosphere  —  the  troposphere.  The  limitations  imposed  by 
these  assumptions  are  flexible  and  depend,  for  instance,  upon  w'hether  the  transmission 
path  under  consideration  is  long  and  near  the  earth’s  surface  or  short  and  liigh  in  the  at¬ 
mosphere.  In  any  case,  the  above  restrictions  are  commonly  accepted  as  being  appropriate 
for  the  problem  of  tropspheric  microwave  propagation  and  lead  to  no  significant  errors  in 
electromagnetic  field  calculations. 

The  methods  described  in  this  report  form  the  basis  of  the  VTRPE  (variable  terrain  radio 
parabolic  equation)  computer  program.^®  The  primary’-  use  of  the  VTRPE  model  is  in  the 
analysis  and  performance  prediction  of  radar  and  communication  systems  operating  in  the 
microwave  region.  In  practice,  the  actual  performance  of  such  systems  in  the  atmosphere  is 
often  quite  different  from  the  characteristics  predicted  based  upon  free-space  propagation. 
Free-space  ranges  are  often  several  orders  of  magnitude  different  from  observed  detection 
ranges  in  the  atmosphere.  The  reason  for  this  discrepancy  is  threefold: 

1.  The  earth’s  surface  is  a  finite  conductor  and  scatters  (reflects)  incident  energy  in 
various  directions,  leading  to  complicated  spatial  interference  patterns. 

2.  The  curved  earth  casts  a  shadow'  giving  rise  to  diffraction  phenomena. 

3.  Inhomogeneities  in  the  atmospheric  index  of  refraction  cause  significtmt  refraction 
or  bending  of  radio  wave  energy. 

To  adequately  represent  the  above  types  of  “anomalous”  propagation  effects,  a  propagation 
model  must  incorporate  full-wave  propagation  physics  and  range-dependent  environmental 
inputs.  The  electromagnetic  paraboDc  wave  eqtiation  does  just  this. 

Tne  remainder  of  this  report  is  divided  into  tw'o  major  sections.  First,  the  fundamental 
physics  and  mathematical  relations  describing  tropospheric  microwave  propagation  over  a 
spherical  earth  are  developed.  Starting  with  the  ba.sic  Maxwell’s  equations  that  govern  the 
electromagnetic  fields,  a  parabolic  wave  equation  is  derived  for  the  nonzero  field  compo¬ 
nents.  Second,  this  parabolic  wave  equation  is  mmierically  solved  by  using  the  split-step 
Fourier  PE  algorithm,  and  the  resulting  v.rror  terms  are  discussed. 
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§2.0  PROPAGATION  PHYSICS 


The  scope  of  this  section  is  to  develop,  from  the  fundamental  equations  of  electrody- 
ncunics,  a  parabolic  wave  equation  (PE)  that  governs  the  propagation  of  electromagnetic 
waves  in  the  troposphere.  In  later  sections,  this  parabolic  wave  equation  wiU  be  solved 
numerically  for  the  electric  and  magnetic  radiation  fields  by  means  of  the  split-step  Fourier 
PE  (SSFPE)  algorithm.  The  model  problem  in  mind  is  the  propagation  of  linearly  polar¬ 
ized  microwave  radiation  through  an  inhomogeneous  troposphere  over  a  spherical  earth. 
The  analytical  development  will  proceed  as  follows:  First,  Max-well’s  equations  and  some 
basic  concepts  of  classical  electromagnetism  are  reviewed  in  order  to  derive  vector  wave 
equations  for  the  electromagnetic  radiation  fields.  By  design,  the  exposition  is  terse,  and 
the  reader  is  referred  to  standard  texts  on  classical  electromagnetic  theory  for  the  de- 
tails.^^"^^  Next,  from  the  vector  wave  equations  an  equivalent  scalar  Helmholtz  equation 
is  derived  for  the  nonzero  field  components.  This  Helmholtz  equation  is  then  transformed 
from  spherical  coordinates  to  an  equivalent  “earth-flattened”  coordinate  system.  Finally, 
the  elliptic  Helmholtz  equation  is  approximated  by  the  Leontovich-Fock  parabolic  wave 
equation. 

The  discussion  is  restricted  to  the  propagation  of  electromagnetic  radiation  in  the  tro¬ 
posphere,  which  is  assumed  to  be  a  linear,  isotropic  nonionized  mediimr.  The  electrical 
properties  of  the  troposphere  are  modeled  by  a  lossy  dielectric.  The  surface  of  the  earth 
will  be  modeled  as  a  nonferrous  dielectric  sphere  with  finite  conductivity.  A  source  of 
electromagnetic  radiation  (i.e.,  a  transmitter)  is  assumed  to  be  located  at  =  (ro,0,0) 
on  the  polar  axis  of  an  earth-based  spherical  coordinate  system  (r,  9,  <f>)  .  The  source  is  as¬ 
sumed  to  emit  linearly  polarized,  monochromatic  radiation  having  a  simple  harmonic  time 
dependence  t  given  by  exp(— ewf),  where  u;  —  2:ir/  is  the  radian  frequency.  (No  generalitji' 
is  lost  by  this  last  restriction,  since  by  Fourier’s  theorem  any  linear  field  of  arbitrary  time 
dependence  can  be  synthesized  from  a  knowledge  of  its  spectral  components.)  Further¬ 
more,  it  is  eilso  assumed  that  the  frequencies  of  interest  are  all  well  above  the  ionospheric 
plasma  frequency  and  that  effects  of  the  earth’s  magnetic  field  can  be  ignored. 


2.1  Maxwell’s  Equations 

With  the  above  restrictions,  the  source-free  monochromatic  Maxwell’s  equations  in  ra¬ 
tionalized  mks  units  are  (’oold  faced  symbols  denote  vector  fields) 


V  X  E(r,u;)  —  -j-?u;B(r,a;), 

V  X  H(r,u;)  =  — iu;D(r,Lt'), 
V-B(r,u;)  =  0, 
V-D(r,a;)  =  0, 


(1) 
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where 


E(r,u;)  =  electric  field  intensit}’'  vector  (V/m), 

H(r  )  —  magnetic  field  intensity  vector  (A/m), 

B(r,  u?)  =  magnetic  field  induction  vector  (Wb/m), 

D(r,aj)  ~  electric  field  displacement  vector  (C/m), 

and  the  impBcit  time-dependence  exp(—iu)t)  has  been  suppressed.  Henceforth,  for  nota- 
tional  simplicity,  the  explicit  frequency  dependence  of  the  field  quantities  is  suppressed. 

The  macroscopic  fields  D  and  H  are  related,  respectively,  to  E  and  B  by  constitutive 
relations  that  characterize  the  electromagnetic  properties  of  the  material  medium  involved. 
For  isotropic  materials  the  constitutive  relations  are 

D(r)  -  £{r)E(r), 

B(r)  =  fi(r)H(r), 

where  e  is  the  absolute  (complex)  permittivity  and  u  is  the  absolute  magnetic  permeability. 
For  most  nonferrous  materials,  the  permeability  /x  is  very  close  to  that  of  fi’ee  space 
Accordingly,  the  medium  permeability  will  be  set  equal  to  its  vacuimi  value:  p(r)  /xq  -  1^ 
also  proves  useful  to  define  the  the  vacuum  wave  number  ko  and  the  dimensionless  relative 
(complex)  permittivity  £  by 

£(r)  = - =£1-1- 1£2  =  1 

fo  ‘*^€0 

A-0  s  uj^^o^o, 

where  eo  is  the  vacuum  dielectric  constant,  a  is  the  medium  conductivity  and  £i  is  the 
usual  relative  permittivity  of  the  medium. 

Maxwell’s  equations  (1)  then  assume  the  form; 


V  X  E(r)  =  4-t'J/xoH(r), 

(2) 

V  X  H(r)  =  — zu;eo^(>')B(r), 

(3) 

/xoV-H(r)-0, 

(4) 

:oV-£(r)E(r)=0. 

(5) 

2.2  Vector  Wave  Equations 

Maxwell’s  equations,  Eqs.  (2)-(5),  8ire  a  set  of  first  order  partial  differential  equations  in 
which  the  E  and  H  fields  are  coupled.  An  equivalent  set  of  second-order  uncoupled  vector 
wave  equations  may  be  derived  from  them,  in  which  the  electric  and  magnetic  fields  occur 
separately.  For  example,  by  taking  the  curl  of  Eq.  (3),  using  standard  vector  identities, 
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and  employing  Eq.  (2),  one  can  eliminate  the  electric  field  and  obtain  an  equation  for  the 
magnetic  field  H  alone: 


V  X  V  X  H(r)  =  — iweoV  x  (£(r)E(r)] , 

=  ^xVxH(r)+(:|£(r)H(r). 


(6) 


By  using  the  vector  identity  V  x  V  x  C  —  V(  V  ■  C)  —  V^C,  and  the  fact  that  the  magnetic 
field  is  solenoidal,  Eq.  (6)  may  be  cast  into  the  form  of  a  vector  wavs  equation: 

V^H(r)  +  X  V  X  H(r)  +  Jl'^e(r)H(r)  =  0.  (7) 

Once  the  magnetic  field  is  found,  by  solving  Eq.  (7),  the  electric  field  is  then  obtained  from 
the  Maxwell  curl  relation  of  Eq.  (3): 


E(r)  = 


i 

Ljeoeir) 


V  X  H(r). 


In  a  similar  fashion,  the  magnetic  field  H  may  be  eliminated  by  taking  the  curl  of  Eq.  (2) 
and  using  Eq.  (3)  to  get  a  vector  wave  equation  involving  the  electric  field  alone: 


V^ECr)  -f-  V 


E(r) 


V£(r) 


Kr)  J 


+  ^'5£(r)E(r)  =  0. 


(8) 


The  Maxwell  curl  relation  of  Eq.  (2)  may  then  be  used  to  obtain  the  H-field: 


—t 


H(r)  = - V  X  E(r). 

wpo 


The  vector  wave  equations  (7)  and  (8)  are  no  more  tractable  than  Maxwell’s  equations 
unless  an  assumption  is  made  regarding  the  field  polarization.  In  many  cases,  the  transmit¬ 
ter  emits  radiation  that  is  linearly  polarized  in  the  meridian  plane  containing  the  receiver. 
For  a  linearly  polarized  transmitter,  the  resulting  electric  field  vector  may  be  assumed 
to  have  components  Ijdng  either  wholly  within  (vertical  polarization)  or  perpendicular  to 
(horizontal  polarization)  the  meridian  plane  containing  the  source  and  observation  point. 
Thus,  for  the  vertically  polarized  case,  the  nonzero  magnetic  field  component  is  H  = 
which  may  be  computed  from  Eq.  (7),  while  for  the  horizontally  polarized  case,  the  only 
nonzero  electric  field  component  is  E  =  which  may  be  obtained  from  Eq.  (8). 


2.3  Vertical  Polarization 

If  the  source  emits  vertically  polarized  radiation,  then  the  electric  field  vector  will  have 
components  in  the  meridian  plane  containing  the  source  and  receiver,  E  —  ErCr  Egig, 
while  the  magnetic  field  vector  has  a  single  nonzero  component  perpenJicvilar  to  this  plane; 


H  =  It  thus  proves  simpler  to  solve  for  the  single  magnetic  field  component  by  using 

Eq.  (7)  and  then  compute  the  remaining  elt?ctric  field  components  from  the  Maxwell  curl 
relation  of  Eq.  (3); 


E(r)  = 
Er  = 
Ee^ 


X  H  =  Er(T)tr  +  EB{v)te, 

i  d{H^sm6) 
tjjeQsr  sin  6  dd 
-i  djrH^) 

iiiCQer  dr 


Working  in  a  spherical  coordinate  system  with  respective  unit  vectors  (cr-  eg.  e^), 

the  first  two  terms  in  Eq.  (7)  become^^ 


,  1  d  ,  . 

r  dr^  r'^smBdB  09 

1  ,  2  coi  e  OH ^ 

^  sin  9  d4>^  sin^  9  sin  9  00 


-2  OH, 


<P 


r^sin^  0^ 


(ir  +  eg  cot  9)  -r 


(9) 


and 


V£ 

£ 


X  VxH  = 


1  Os 


e  d<l>  sin  9 


e^  1 
r  [  r  sin  9 


d{rH^)  ,  1  d{sm9H<f,y 

^  Or  ^  sin  ^  09 

djsmBH^)  1  Oc  IdeOirH^y 
09  £  09^  £  dr  dr 


(10) 


Equations  (9)  and  (10)  reveal  that,  in  general,  the  indi\'idual  vector  components  of  the 
magnetic  field  vdll  be  coupled  due  to  the  presence  of  the  ^  term.  This  coupling,  or  depo¬ 
larization,  is  analogous  to  the  Faraday  rotation  observ'ed  in  optically  active  materials.  A 
significant  simplification  arises  if  the  propagation  meditun  is  a.ssunied  to  have  no  compo¬ 
nents  of  Ve(r)  parallel  to  the  source  magnetic  field  (i.e.,  •  Vs  =  0).  In  this  case,  Eqs.  (9) 

and  (10)  will  have  only  a  single  nonzero  component: 


1  d^{rHft>) 


4  •  —  X  V  X  H  =  ~ 


+ 


^(sin.^^^ 


) 


Hr 


.2  c;r,2  0 


sin 9  do''  09  r*  sin 

1  dism9Ho)ld£  1  Os  d{rH^) 

r  sin  9  09  e  09  e  Or  Or 


(11) 

(12) 


The  r  and  ^-deri'.’atives  appearing  in  Eqs,  (11)  and  (12)  may  be  combined  to  give  the 
followng  compact  form  for  the  vector  wave  equation  (7): 


eOlOirH^)  e  0  smBOJI^  ,  f ,2  1  cot9d£\  _ 

r  Or  e  dr  r‘^  sin  9  06  e  09  \  sin“  #  r-e  OBJ 
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As  it  stands,  Eq.  ( 13)  is  not  in  the  desired  Helmholtz  form  but  can  be  transformed  into 
it  by  a  change  of  the  dependent  variable; 


n(r)  =  VeOr)- 


(14) 


where  n  is  the  complex  index  of  re&’action  referenced  to  the  vacuum  dielectric  constant  eo- 
In  terms  of  u,  the  partial  derivative  terms  in  Eq.  (13)  can  be  shown  to  have  the  form 


^  1  d(rH^)  ^  I  d  I  d(u\/Te) 
dr  c  dr  dr  e  dr 


I  d  du 
r  dr  dr 


r  1  d  du 
I  r  dt' 


and 


d  sin9  dH^  _  ^  ^  sin^  d{u\/€ esc 9} 

dd  e  d6  y/rsmOdQ  e  dd 


/sin^ 

d^u  j 

( 1 

csc^  6 

cot  6 

de 

I  d‘^£  2  \ 

=  V  re 

r 

4  + 

2e 

de^ 

(15) 


(16) 


Now,  with  the  new  field  vEiriable  Ur,  Eq-  (13)  takes  on  the  desired  Helmholtz  equation 
form: 


1  d  duv  1  d^Uv 
r  dr  dr  ^  d9^ 


+  [A'0‘(*’)  +^^t'(r)]wr 


=  0, 


(17) 


with 


3  csc^  ^  cot  Ode  id^e~2  ^2  gf^s~2 

4  2r^£  d9  ^  dr^  r^  dO^ 
3csc^(?  cot^^  n  d^n~^  d^n~^ 

4  r^n  dd  r^  d6^  dr^ 


In  Eq.  (17)  the  full  d^dependence  is  retained  in  the  e-dependent  terms  in  the  coefficient 
of  Ut.,  even  though  ^  =  0  in  the  deri%'ation.  This  appears  to  be  a  contradiction,  wliich 
will  now  be  justified,  albeit  heuristically. 

The  terms  involving  ^  occur  in  the  coefficients  of  the  partial  deri\'ative  operators  in 
Eq.  (11)  and  will  be  neglected.  Basically  this  amounts  to  neglecting  small  variations  in  the 
field  amplitude,  but  retaining  the  <j>-dependence  in  the  term  which  affects  the  phase 
of  the  field.  This  is  most  easily  understood  in  the  context  of  geometrical  optics,  where  the 
large  /tq  solution  to 

M  -f-  kle(r) 


0, 


is  ex|)rea5ecl  in  tile  form^^ 


H4,  = 

where  the  phase  function  S  satisfies  the  eikonal  equation 

(VSf  =  £(r). 

Thus,  small  changes  in  s  can  lead  to  large  overall  phase  errors  in  the  eikonal  k^S  and  cause 
appreciable  errors  in  the  field.  This  is  not  the  case  for  the  amplitude  function  Q,  which  is 
only  weakly  dependent  on  variations  in  e. 

Neglecting  the  azimuthal  variation  in  Eki.  (10)  but  including  it  in  the  k^e  factor  is 
termed  a  2|-D  approximation  to  the  fiili  3-D  Helmholtz  equation.  It  heis  the  advantage  of 
retaining,  to  lowest  order,  the  effects  of  azimuthal  enviroimieiital  v'ariability  wliile  reducing 
computational  complexity. 

2.4  Horizontal  Polarization 


If  the  source  emits  horizontally  polarized  radiation,  then  the  electric  field  vector  hzis  the 
single  non25ero  component  E  =  Analogous  with  the  case  of  vertical  polarization, 

the  single  electric  field  component  may  be  obtained  from  Eq.  (8)  which,  in  spherical 
coordinates,  becomes 


ia2(rE,(r))  Id.  ,S£*(r)  , 


Again  paralleling  the  vertical  polarization  case,  Eq.  (19)  can  be  converted  into  a  scalar 
Helmholtz  equation  by  introducing  a  new  dependent  field  variable  u/,,  defined  by 

u/i  =  «ft(r)  =  \/rsin^£’^(r),  (20) 

and  dropping  the  (;^-derivative  term  in  Eq.  (19).  With  the  new  field  variable,  Eq.  (19)  then 
becomes 

1  d  duh  1  d^uh  r,2  ,  ,  r  /  \1  r 

^  1^'W  +  nh  -  0,  (..1) 


where 


i>eh{T,B) 


4r‘^  sin^  9 


Thus  both  the  horizontal  ^lnd  vertical  polarization  cases  can  be  expressed  in  terms  of  a 
scalar  Helmholtz  equation  of  the  form 

+fc(r)j  u(r)  =  0,  (23) 

where  u  and  6c  ai-e  defined  by  Eqs.  (14)  end  (18)  for  vertical  polarization,  and  by  Eqs.  (20) 
and  (22)  for  horizontal  polarization. 
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The  additional  term  Se  is  generally  small,  except  in  two  cases;  The  first  occurs  for 
vertical  polarization  with  propagation  near  the  boundaiy^  separating  dissimilar  dielectrics 
with  finite  conductivity,  in  which  case  the  spatial  derivatives  of  the  index  of  refraction  n 
may  become  large.  The  second  occurs  for  long  propagation  paths  at  low  frequencies  where 
the  wave  number  kQ  is  small  and  then  the  spherical  correction  terms  become  important. 
Inclusion  of  the  6s  correction  allows  one  to  deal  in  a  unified  manner  with  propagation  over 
idealized  bormdaries  as  well  as  treat  propagation  in  the  presence  of  subsurface  overbmden. 

2.6  Earth-Flattening  Transform 

The  coupling  of  r  and  6  variables  in  Eq.  (23)  is  not  desirable  for  numerical  purposes 
and  can  be  eliminated  by  means  of  a  simple  transformation  from  the  spherical  (r, 
coordinate  system  to  an  equivalent  cartesian  cooiuinate  system  defined  by 

C  =:aln(r/a), 

(~aBcos(f>,  (24) 

Tf  =  aO  sin 

where  a  is  the  earth’s  radius,  and  C  =  0  corresponds  to  the  surface  of  the  earth.  This 
transformation,  suggested  by  Pryce,^®  is  known  as  an  earth-flattening  transformation,  and 
with  the  approximation  ss  ft  was  also  studied  by  Pekeris.^^ 

In  terms  of  the  new  coordinates  defined  by  Eq.  (24),  the  metric 

ds^  =  dr^  -f  r^dB^  -f  sin^ 


becomes 


ds2  = 


d^  +  drf'  +  dQ^  —  (1  — 


6>2  ’  + 


(25) 


In  most  problems  of  interest,  ^  >C  1,  in  v/hich  c^lse  Eq.  (25)  can  be  expanded  in  powers  of 
B  to  yield 


d.'?  = 


d^  -|-  dr}^  +  — 


2  iidr}  -  rjdO' 


3a2 


+  0(a-") 


(26) 


It  is  a  fairly  simple  matter  to  determine  the  geodesics  of  the  non-Euclidean  (i.e.,  non¬ 
fiat)  space  whose  metric  is  given  by  Eq.  (26).  Let  the  source  point  be  located  at  (^  =  0, 7  = 
0,  ^  =  ^0)  limit  attention  to  geodesics  w'hich  pass  through  the  source.  From  symmetry', 
it  is  obvious  that  the  geodesics  lie  in  a  plane  containing  the  We  may  without 

any  loss  of  generality  confine  ourselves  to  geodesics  in  the  ^(^-plane.  These  geodesics  will 
then  coincide  with  straight  lines  through  the  source,  the  equations  of  which  are  given 
parametrically  by 


^  =  a  tan  ^  ^tanV’  +  —  sect(ie  —  arp, 

C  =  Co  +  ^aln  fl  -{-  —  sind’e“‘’°/®  -t-  , 

2  \  a  / 


(27) 
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where  s  is  the  geodesic  distance  (arc  length)  from  the  source.  The  parameter  ip  is  a,  constant 
along  the  geodesic,  being  the  angle  the  geodesic  make.s  with  the  horizontal  at  the  source. 
In  terms  of  the  original  spherical  coordinates  <p),  Eq.  (27)  reduces  to  the  equation  for 
a  straight  line; 


r  sin  0  =  s  cos  ip, 
r  cos  &  =  ssinip  +  . 


Also,  the  distance  along  the  geodesic  between  the  points  (0,0,  Co)  and  (^,0,Q  is  simply 

1/2 


s  =  a 


g2<o/a  _  2e(<«+^)/“cos(C/a)  + 


In  what  follows,  the  meridian  plane  containing  the  source/receiver  will  be  taken  to  be 
the  r 6-plane.  In  this  case,  the  earth-flattening  transformation,  Eq.  (24),  is  simply  defined 
in  terms  of  new  “cartesian”  coordinates  {x,z)  by 


X  =  a6, 

h 

z  ~  aln(l  +  h/a)  Rs  h(l  —  — ),  if  h/a  C  1 

^iX 


(28) 


where  h  —  r  ~  a  is  the  local  altitude  referenced  to  the  mean  earth  radius  a.  This  leads  to 
the  following  Helmholtz  separated  form  for  Eq.  (23)  : 


■  0^ 
dx^  ^  dz^ 


+  K\x,z) 


w(x,  z)  =  0. 


(29) 


The  new  dependent  N’ariable  w  and  the  “effective”  wave  number  K  are  defined  for  horizontal 
polarization  by 

wfi  =  VrsiiiBE^ir), 

e‘*‘^/^*'^“^y'osin(r/a)E0(x),  (30) 

Ks  y/x  if  x/a  <C  1  and  z/a  <C  1 


and 


KJ, 


1.2  ._2 
K(\ 


071 


3csc^(a'/o) 
4^2  ■ 


if  x/a  <?:  1 


with  m  the  modified  index  of  refraction,  defined  by 


(31) 


m  =  n(r,6,  (^)-  =  n(r)(l  +  -). 

a  a 


(32) 
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For  vertical  polarization,  the  new  variables  are 


and 


\/rsm$ 

Wv  = - HJr), 

n 


^^+3z/(2.)\Aplp„ 

m(ar,  z)  ^ 


m 


if  x/a  <C  1  and  zja  <C  1 


K(;{x,z)  =  - 


2  2  3csc^(ar/a)  cot{x/a)dm 


4a^  am 

d^m~^  mdm~^ 


—  7n- 


5^2 


m- 


5z2 


a  ^2 


(33) 


(34) 


If  the  medium  is  spherically  stratified  and  has  no  range  dependence,  then  the  effective 
wave  number  takes  the  form 


A’2(r,  z)  klmlffiz)  - 

1 

"^e//(^)  =  rn^{z)->r  ~p^ 

«  u^(2)(l  +  22/a),  if  2/a  <  1 


1  d‘^m  2 

m  ^22 


if  r/a  -C  1 
2 

a. 

am  dz 


1  dm\^  1  dm 
dz  ) 


which  is  the  same  form  as  the  earth-flattening  approximation  discussed  by  Pekeris.  The 
Pekeris  transform,  however,  was  based  upon  2  =  r  —  a  and  is  only  valid  for  small  2/a  in 
contrast  to  that  given  by  Eq.  (28),  which  is  exact. 


2.6  Boundary  Conditions 

The  solutions  to  Maxwell’s  equations  are  not  unique  until  boundary  conditions  on  the 
fields  are  prescribed.  These  boundary  conditions  are  (1)  a  Sommerfeld  radiation- type 
boundary  condition  at  infinity^® 


^h^r  -  ik^Aj  ,  (35) 

where  A  denotes  the  field  component  (i.e.,  A  =  for  horizontal  polarization,  >1  =  for 
vertical  polarization),  and  (2)  continuity  of  the  tangential  electric  and  magnetic  fields  at 
the  earth’s  surface,  r  =  a.  If  the  propagation  medium  is  lossy  —  i.e.,  the  dielectric  constant 
is  complex  with  a  positive  imaginarj'  part  ~  then  the  Sommerfeld  radiation  condition  may 
be  replaced  by  the  simpler  condition  that  the  fields  vanish  as  r  00  . 

The  continuity  of  tangential  field  components  at  r  =  a  is  implemented  by  modeling  the 
earth  as  a  locally  homogeneous  dielectric  with  finite  conductivity.  This  is  a  reasonable 
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approximation  at  microwave  frequencies,  since  the  penetration  depth  of  the  fields  into  the 
earth’s  surface  is  small  compared  to  spatial  variations  in  the  surface  dielectric  properties. 

For  example,  in  the  case  of  microwave  propagation  over  seawater,  it  can  be  shown  by 
standard  techniques  that  the  electric  and  magnetic  fields  within  the  ocean  decay  exponen¬ 
tially  with  distance  from  the  air-sea  boimdary.^^  The  vertical  scale  length  over  which  the 
field  components  in  the  water  decay  to  1/c  of  their  value  at  the  interface  is  termed  the 
skin  depth  6.  The  high-frequency  form  of  6  is^^ 

S  ~  y/2/fiocrsU}, 


where  <Ta  is  the  surface  conductivity.  A  typical  value  for  the  conductivity  of  seawater 
is  <T  «  4  S/ni,  which  means  the  order  of  magnitude  of  S  is  about  250/ \/f,  with  /,  the 
frequency,  in  hertz.  For  typical  marine  radar  frequencies  (/  ss  10^®  Hz),  S  k  0.025  m  . 

The  small  skin  depth  means  that,  as  an  approximation,  the  boundary  condition  of 
tangential  field  continuity  can  be  replaced  by  the  Leontovich  surface  impedance  boundary 
condition^® 


d{rA) 


dr 


=  -Z{rA) 


r—a 


r~a 


(36) 


The  quantity  Z  is  the  local,  flat  surface  impedance  given  by  Fresnel  formulas: 


ikQ 

es 


cos 


for  horizontal  polarization, 
for  vertical  polaiizatioa, 


with  £3  the  complex  dielectric  constant  of  the  earth’s  surface  and  the  grazing  angle 
of  the  field  at  the  earth’s  surface.  In  most  applications,  the  surface  grazing  angle  is  verj’^ 
small,  and  it  is  appropriate  to  use  the  limiting  forms 


Z/* (37) 
Zt,  =  — (38) 

^3 


The  Leontovich  impedance  boundary  conditions,  Eq.  (37)  and  Eq.  (38),  may  also  be 
transformed  to  earth-flattened  cartesian  coordinates  where  they  become 

i(;/,(x,0)  (39) 

for  horizontal  polarization  and 

-fR'O— - ^Vi)„(r,0)  (40) 

s=0  / 

for  vertical  polarization, 


du\,{x,  z) 


dz 


1  dm{x^z) 
2a  m(r,0) 


dz 


dwhix,z) 

dz 
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2.7  Propagation  Factor 


In  analyzing  the  'rarious  propagation  phenomena,  it  is  useful  to  separate  those  system 
parameters  not  influenced  by  the  environment,  antenna  radiation  characteristics  for  exam¬ 
ple,  from  propagation  effects  that  are  environmentally  influenced,  such  as  ducting  caused 
by  spatial  variations  in  the  refractivity.  Following  Kerr,^^  define  the  one-way  generalized 
transmission  equation,  wliich  relates  the  power  received  by  an  omnidirectional  receiver  at 
a  point  in  space,  Pr(r),  to  the  power  emitted  from  a  transmitting  antenna,  Pt,  by 


Pr{r) 

Pt 


Gt 


12 


L(2A:o/?)J 


(41) 


where 


Gt 

ko 


—  transmitting  antenna  power  gain, 
=  vacuum  wavenumber  = 


R  =  distance  between  transmitter  and  receiving  point, 
F  =  pattern  propagation  factor. 


The  radiating  characteristics  of  an  anteima  are  specified  in  terms  of  the  antenna  radia¬ 
tion  pattern  function  f{0,  <f>)  where  {$,  <f>)  are  the  zenith  and  azimuthal  angles,  respectively, 
of  a  spherical  coordinate  system  centered  at  the  antenna,  with  polar  axis  pointed  in  the 
direction  of  maximum  transmission.  The  antenna  radiation  pattern  function  /  is  defined 
to  be  the  ratio  of  electric  (or  magnetic)  field  strength  E(^,  4>)  radiated  in  the  direction 
(^,  0),  to  the  peak  transmitted  field  strength  Eq: 


E(g,0) 

Eq 


The  radiation  pattern  function  is  related  to  the  time-averaged  Poynting  vector  S  of  the 
radiated  wave  field  by 

S[eA)  =  \f{eA)?SQ^ 

where  Sq  is  the  energy  flow  per  unit  area  corresponding  to  the  peak  field  Eq  .  In  general, 
the  antenna  radiation  pattern  function  /  is  a  complex  valued  quantity. 

The  antenna  gain  G  is  expressed  in  terms  of  the  pattern  function  as 


G  = 


Note  that  the  radiation  pattern  function  /  and  the  corresponding  anteima  geiin  G  are 
defined  with  respect  to  free-space  propagation  conditions,  and  therefore  do  not  include 
any  environmental  effects. 
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Environmental  effects  are  included  in  the  generalized  transmission  equation  via  the 
pattern  propagation  factor  F,  The  pattern  propagation  factor  is  defined  as  the  ratio  of 
the  field  magnitude  at  a  point  in  space,  E{r),  to  the  magnitude  of  the  field  at  the  same 
point  but  under  free-space  conditions,  EoC*")- 


E(r) 

Eo(r) 


(42) 


This  definition  of  the  pattern  propagation  factor  assumes  that  the  transmitting  antenna 
is  aligned  with  its  maximum  response  axis  pointed  directly  at  the  observation  point.  The 
pattern  propagation  factor  is  the  key  part  of  Eq.  (41),  and  is  the  fundamental  quantity  to 
be  computed  by  the  PE  model.  Because  of  the  large  dynamic  ranges  of  the  various  terms 
in  Eq.  (41),  it  is  customarj'^  to  work  in  dB-space  by  defining  the  dimensionless  variables 
PF  and  PL  by 


PF  =  20log|El,  (43) 

PL  -  201og(2feoi?)  -  20  log  IFl.  (44) 

The  quantity  PF  is  often  called  the  propagation  factor,  w'hile  PL  is  denoted  the  path  loss. 

The  remaining  task  is  to  compute  the  pattern  propagation  factor  F  defined  by  Eq.  (42), 
which,  for  linearly  polarized  radiation,  is  just  the  ratio  of  the  <;&- components  of  the  electric 
or  magnetic  field  vectors; 


F^FUr)- 

horizontal  polarization 

(45) 

vertical  polarization 

(46) 

where  and  are  the  free-space  electric  and  magnetic  field  components.  By  con¬ 
vention,  the  pattern  propagation  factor  is  defined  with  respect  to  the  free-space  field  of  a 
unit-strength  point  dipole. 

Following  Papas,^^  the  free-space  dipole  fields  are  computed  by  using  the  dyadic  Green’s 
function  r(r,  Tq),  which  is  the  solution  of 

V  X  V  X  r(r,ro)  -  fc^r(r,ro)  =  u6{r  -  fq), 

where  u  is  the  unit  dyadic.  The  dyadic  Green’s  function  F  can  be  expressed  in  terms  of  a 
scalar  Green’s  function  Gq  as 

r-(u4-pVV)Go, 

where  the  scalar  Green’s  function  Go  satisfies 


(V^-hfc2)Go(r,ro)  =  -6(r-ro). 


(47) 
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The  appropriate  outgoing  wave  solution  of  Eq.  (47)  is 


Go(r,ro) 


1 

47r  |r-ro| 


(48) 


Two  types  of  dipole  fields  are  required:  (1)  a  vertically  polarized  field  arising  from  a 
point  vertical  electric  dipole  (VED)  with  electric  dipole  moment  oriented  along  the  s-axis, 
p  =pgig;  and  (2)  a  horizontally  polarized  field  arising  from  a  point  vertical  magnetic  dipole 
(VMD)  with  magnetic  dipole  moment  oriented  along  the  2-axis,  m  =  rriztz.  Employing 
the  dyadic  Green’s  function  F,  the  electric  and  magnetic  fields  from  a  VED  are  given  by 

E,j(r)--ivx(px  VGo),  (49) 

Hed(r)  =  iujp  X  VGo,  (50) 

while  the  fields  from  a  VMD  are  given  by 

Hmd(*-)  = -V  X  (m  X  VGo),  (51) 

Em<f(r)  =  x  VGq.  (52) 


For  both  the  VED  and  VMD  cases,  let  a  unit-strength  point  dipole  be  located  at  ro  = 
(ro,0,0)  in  a  spherical  coordinate  system  with  the  dipole  moment  oriented  along 

the  polar  axis.  Using  the  scalar  free-space  Green’s  function  Gq,  Eq.  (48),  and  equations 
(50)  and  (52),  the  (^-components  of  the  dipole  fields  are  given  by 

«e,(r)  =  (^)  (to  +  i)  .  VED  (53) 

(^)  (to  +  i)  ,  VMD  (54) 


where  i?  =  jr  —  rol  =  y  4  —  2rro  cos  6,  and  6  is  the  polar  angle. 

The  propagation  factor  for  horizontal  polarization  F  =  F/,  then  becomes 


Fh{r)  = 


E^{r) 


Emdir) 
47r 


u(r)li?^ 


(rsin^)2 

while  that  for  vertically  polarized  radiation  F  =  Ft,  is 


1  4  ihR) 


-2 


Ev{r) 


E<p{r) 

’ 

47r  jn(r)a(r)|i?^ 
fepw  (rsin^)2 


1  -t-  ikoR)~^ 


1 

2 


(55) 


(56) 
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To  express  the  propagation  factor  F  in  earth-flattened  (r,  z)  coordinates,  the  dipole 
fields  defined  in  Eq.  (53)  and  Eq.  (54)  need  to  be  converted.  First,  the  source- receiver 
separation  distance  R  becomes 

=  Ir  -  rol  =  yj (r  -  r©)^  -f  4rT*o  sin2(^/2), 

=  Jsinli2(l^)  +  sm2(^), 

V  2g  2a 

«  (1  +  -  zaf  +  0(a-^/a^), 

-t  yj{z  -  -f  if  xja  <C  1  and  {z  -f-  2o)/«  ■C  1  .  (57) 

The  last  relation  in  Eq.  (57)  has  a  relative  error  of  <  2.5x10”^  for  horizontal  ranges  x  <  200 
km  and  altitudes  2  <  30  km.  In  like  manner,  the  angular  dipole  coefficient  (r/J?,)sin^  has 
the  limiting  forms: 


•—  sin  ^  =  {a/R)e^^^  s\n(x/a). 

R 

X 

—*■  -=,  if  xfa  -C  1  and  (2  -|~  2o)/o  ^  1 
R 

1,  if  (2  -  zo)/x  <  1  ,  (58) 


Finally,  the  propagation  factor  takes  the  following  limiting  form  for  altitudes  and  range's 
typical  of  tropospheric  propagation: 


and 


F/,(x) « 


4ffi?2  _3  |tr’(x)l 

_ 2  *  ^  - 

y/T+ThW^' 


4ity/x  |«’(^)l 

v^l  -f  (fcor)-^  ’ 


if  (2  -  2o)/x  <  1 


Fv{x) 


4t:R^  _3  |7n(x)u»(x)l 

— - — X  2  *  -  ~ 

^0^'  y/i  -^{koR)™^ 


47r  iyx|m(x)tt'(x)| 

-^1  -I-  {kQx)~^ 


if  (2  —  zo)/x  <C  1  . 


(59) 


(60) 


2.8  Parabolic  Wave  Equation 


Although  Eq.  (29)  follows  exactly  from  Maxweirs  equations,  analytic  solutions  are  only 
possible  if  K  is  independent  of  x.  Nor  is  direct  numerical  solution  an  easy  alternative. 
This  follows  from  the  fact  that  Eq.  (29)  is  an  elhptic  partial  differential  equation;  hence 
solving  for  a  solution  at  one  point  requires  that  it  must  be  numerically  solved  over  the 
entire  propagation  domain  simultaneously.  Standard  finite  difference  approaches  to  solving 
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elliptic  boundary  value  problems  also  tend  to  require  griding  sdiemes  on  the  order  of  a 
wavelength.  For  the  long  ranges  involved  in  tropospheric  propagation,  this  would  lead  to 
unacceptably  large  problem  sizes. 

To  deal  with  these  issues,  Leontovich  and  Fock  introduced  the  concept  of  a  parabolic 
wave  equation  to  treat  the  problem  of  transhorizoii  radio  wave  propagation.  Their  approach 
may  be  summarized  as  follows:  First,  an  envelope  transformation  is  performed  to  remove 
the  anticipated  rapid  phase  variation; 


w(x)  = 


(61) 


where  k  is  a.  reference  wave  ntunber.  (  For  most  tropospheric  problems,  the  natural  choice 
is  k  =  ko.)  Next,  Eq.  (61)  is  substituted  into  Eq.  (29)  to  get 


2ik 


dx 


+ 


dz^ 


?/)(x)  = 


dx^  ■ 


(62) 


The  essence  of  the  Leontovich-Fock  PE  approximation  is  to  assume  that  variations  of  the 
envelope  function  with  range  are  small  relative  to  other  terms,  and  therefore  to  drop  the 
term  to  get  a  parabolic  equation  in  the  longitudinal  coordinate  x: 


(63) 


The  advantage  of  Eq.  (63)  over  Eq.  (62)  is  that  it  is  parabolic  rather  than  elliptic  in 
range.  This  means  that  it  is  an  initial  value  problem  and  can  be  solved  efficiently  by 
marching  algorithms.  Also,  the  surface  boundary  conditions  retain  the  same  form  in  the 
PE  approximation  as  in  the  original  Helmholtz  equation. 

The  justification  for  dropping  the  second  deri\'ative  term  is  usually  made  on  the  groimds 
that  the  envelope  fimction  is  a  slowly  varying  fimction  of  the  range  coordinate  x.  While 
pedagogically  correct,  such  a  derivation  of  PE  does  not  yield  any  insight  into  the  approx¬ 
imation  nor  the  errors  that  are  incurred  in  using  it.  An  alternative  approach  to  deriving 
the  PE  is  via  factorization  of  the  elliptic  wave  equation.  If  we  define  an  operator  Q{x)  by 


dz‘^ 


(64) 


then  Eq.  (29)  can  be  expressed  in  the  equivalent  factored  form 


)  +  j 


vi{x)  —  0. 


(65) 


(The  notation  ~  the  commutator  of  the  operators  ^  and  Q.) 
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For  range-independent  propagation,  the  commutator 
satisfied  by  out%vardly  propagating  waves  is  just 


=  0  and  the  equation 


=  iQ(x)w(x).  (66) 

Equation  (66)  is  an  exact  formal  solution  to  the  range-  independent  Helmholtz  equation, 
and  is  the  most  general  PE  that  is  exact  for  range  independent  media.  Following  Tappert,^ 
a  parabolic  wave  equation  with  the  Q  operator  defined  by  Eq.  (64)  will  be  denoted  as  the 
general  PE  (GPE)  and  the  Q  operator  will  be  denoted  as  the  GPE  propagator.  GPE  is  the 
most  complete  PE  that  is  evolutionary  in  range  and  neglects  backscattering.  For  range- 
independent  environments,  it  is  exact  within  the  limits  of  the  far-field  approximation,  and  is 
the  starting  point  for  all  numerical  PE  algorithms.  In  following  sections,  the  computational 
techniques  ased  to  solve  the  PE  will  be  derived. 


2.9  Variable  Terrain  Physics 


The  preceding  sections  have  dealt  with  the  problem  of  radio  wave  propagation  over 
a  smooth  spherical  earth.  This  enabled  the  surface  boundary  conditions  on  the  electro¬ 
magnetic  field  components  to  be  satisfied  on  the  conformal  spherical  shell  r  =  a,  where 
a  is  the  average  earth  radius.  In  reality,  the  surface  of  the  earth  is  not  smooth  but  has 
a  nonunTorm  topography.  What  effect  does  this  varying  terrain  have  on  the  propagation 
of  radio  vv-aves?  Intuitively,  one  expects  that  the  Leontovich  surface  impedanc.  boundary 
condition 


d{rA) 


dr 


=  -2(rA) 


r=a 


r=ra 


(67) 


will  be  modified.  This,  in  fact,  is  the  case. 

If  the  terrain  relief  {^{0)  is  a  fvmction  of  the  meridian  angle  9,  such  that  ."=0  4-  ^{B)  is 
the  radial  distance  corresponding  to  the  local  surface,  the  boundary  condition  in  Eq.  (67) 
will  become 


^(rA) 

5n 


=  -Z(rA){ 

r=a4-<(S)  lr=a-|-C(^) 


(68) 


where  A  denotes  the  electromagnetic  field  (^-component  (i.e.,  A  ^  for  horizontal  polar¬ 
ization;  A  —  for  vertical  polarization),  and  -  is  the  normal  derivative  to  the  surface 
defined  by  This  surface  normal  derivative  is  specified  by 


d 

dn 


a 

57 


1  5 
..  M 


Vi  +  7^  ’  r  dB 


(69) 


Using  Eq.  (69),  the  Leontovich  boundary  condition  may  be  written  as 


d(rA) 


dr 


—  T 


r=a-(-C 


dA 

dB 


~a+( 


Z  v^l  +  r2(rA) 


r=a-|-(; 


(70) 
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The  new  boundarj^  condition,  Eq.  (70),  is  seen  to  be  more  complicated  than  the  smooth 
surface  condition  in  Eq.  (67)  due  to  the  presence  of  the  term  involving  This  com¬ 
plication  is  difficult  to  deal  with  numerically,  and  it  seems  logical  to  investigate  whether 
a  modification  to  the  standard  earth-flattening  transformation  would,  in  elfect,  “smooth” 
out  the  variable  terrain.  Accordingly,  let  us  define  the  surface-locked  cartesian  coordinate 
system  (a:,  ft),  where  x  is  the  local  horizontal  (range)  coordinate  and  ft  is  a  new  verti¬ 
cal  coordinate  measured  with  respect  to  the  local  terrain  elevation.  The  new  flattening 
transform  is  defined  by 

X  —  ad, 

h  ~  aln{r/a)  —  Ci^). 

In  terms  of  the  (a.-,  ft)  coordinates,  the  normal  derivative,  Eq.  (69),  takes  on  the  simple 
form 

A- A 

dn  dh 

and  the  boundary  condition  Elq.  (70)  becomes 


d{rA) 


dh 


h-Q 


=  -Z{rA) 


/i=0 


Thus,  the  new  transformation  defined  by  Eq.  (71)  has  indeed  led  to  a  simpUfication  of  the 
surface  boundarj'^  condition — but  at  what  price? 

If  the  transform  defined  by  Eq.  (71)  is  applied  to  the  scalar  Helmholtz  equation  in 
spherical  coordinates 


Qj.2  ^  r  dr 


+ 


1  52 


4-  kQn^{r,d)  +  Sn{r) 


u{r,e)  =  0, 


(72) 


then  the  result  is 


d^u 
9ft  2 


2a  d^u 
a  dxdh 


ft  du  d^u 


+ 


0  0 

k^m  -f  Sm 


ul  h,  r)  =  0 


(73) 


where 

_  5C  .  _  92c 

“~9^’  ““9^2- 

Equation  (73)  is  no  longer  in  Helmholtz  format — the  flattening  transformation  has  de¬ 
stroyed  its  separabilitj'  by  admixing  the  xft-deri\'atives. 

The  complications  in  Elq,  (73)  can  be  shown  to  arise  from  the  term  in  Eq.  (72). 
Ultimately  the  Leontovich-Fock  parabolic  approximation  is  to  be  applied,  w'hich  amounts 
to  keeping  terms  to  first  order  in  Tliis  suggests  an  approach  that  first  makes  a  PE 
approximation  in  spherical  coordinates  to  Eq.  (72)  and  then  applies  the  earth-flattening 
transformation  Eq.  (71).  Accordingly,  let  us  define  the  envelope  function  by 


u{r,e}  = 
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and  then  write  Eq.  (72)  in  parabolic  form  as 


1  ^  ^  4.  ^ 

r  dr  dr  r  dO 


( 


knTl^  +  671  — 


Next,  let  us  apply  the  transform,  Eq.  (71),  to  get 


d^7l7 


4-  i2k 


\dx 


I 

The  term  involving  can  be  removed  by  a  simple  change  of  variable; 


(74) 


(75) 


7p  =» 

in  which  case  Eq.  (75)  takes  on  the  final  form: 


d^tp 


■  r  dip 

+  - h 

dx 


IrQrn^  +  ~  2hd  —  1) 


ip  —  0. 


(76) 


Comparing  with  the  form  of  the  PE  for  a  smooth  sm-face,  Eq.  (63),  the  effects  of  variable 
terrain  are  seen  to  be  incorporated  into  a  modified  wavenumber  defined  by 


=  ^0^^  +Sm  +  P  2hd  -  . 


(77) 


The  inclusion  of  variable  terrain  in  the  PE  approximation  is  thus  modeled  by  modifying 
the  medium  index  of  refraction  via  Eq.  (77)  and  then  proceeding  with  numerical  solutions 
as  in  the  smooth-surface  PE.  Having  solved  Eq.  (76)  for  the  field  in  xh-coordinates,  it  is  a 
simple  matter  to  transform  them  back  to  physical  space  and  then  compute  the  propagation 
factor  F. 


§3.0  NUMERICAL  SOLUTIONS 


This  section  describes  methods  for  the  numerical  solution  of  the  generalized  one-way 
parabolic  wave  equation 

dihix.) 

=  (78) 

where  the  generalized  PE  (GPE)  operator  Q  is  defined  to  be 


g(x)  = 


yd^ 


+  KHx)  - 


(79) 


with  k  a  reference  wave  number.  Because  of  the  square  root  appearing  in  Eq.  (79),  the 
usual  techniques  for  solving  parabolic  partial  differential  equations  are  not  possible.  This 
is  because  Q  belongs  to  the  class  of  pseudodifferential  operators^^  since  it  contains  both 
multiplicative,  and  differential,  operators  under  the  radical.  Hence,  Q(x)i/’(x) 
cannot  be  expressed  as  a  finite  Taylor  series  in  local  operators  about  the  point  x.  The  Q 
operator  may  be  expressed  in  terms  of  a  dimensionless  '‘kinetic  energy”  operator  T  and  a 
“potential  energ3’”  operator  V  as 


Q{x)  =  frov^l“T(z)-V(r,z)  -  k , 


(SO) 


where  the  kinetic  and  potential  energy  operators  are 


T{z) 


dz^  ' 


I^(x,5)  =  1  -  K^{x,z)/kQ  . 


(81) 


If  j|r  +  V''||  ■<  1,  then  a  formal  Taylor  series  expansion  of  Eq.  (80)  gives 

Q(x)  «  Qiix)  =  kQ-k~  ^{T  +  V),  (82) 

which  is  the  standard  parabolic  equation  (SPE)  operator  first  proposed  by  Tapper!.^  The 
SPE  approximation  to  Q  basically  assumes  that  the  propagation  occurs  within  a  small 
cone  of  angles  centered  about  k. 

An  alternative  approximation  to  Q  was  developed  by  Feit  and  Fleck'^'^  for  propagation  in 
optical  fibers  and  used  by  Thompson  and  Chapman^“  for  imderwater  acoustic  propagation. 
In  this  form,  the  GPE  operator  is  approximated  as: 


Q(x)  Q2{x)  =  A-o(v/r^  -  1)  -1-  A'(x)  -  k .  (83) 

This  approximation  is  known  as  a  wide-angle  parabolic  equation  (WAPE)  operator  since 
it  is  v'alid  for  much  wider  propagation  angles  than  the  SPE. 
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Suppose,  for  the  moment,  that  the  GPE  Q-operator  is  not  a  ftinction  of  the  range 
coordinate  x:  Q  —*  Q{z).  Then,  the  formal  solution  of  Eq.  (78)  is  just 

il'(x,z)  ~  z)  , 

where  the  exponential  of  the  operator  Q  is  defined  by  its  power  series  expansion.  The 
conventional  nmnerical  approach  to  solving  Eq.  (78)  is  to  use  finite-difference  methods  and 
approximate  the  exponential  operator  exp[iQ(a;  —  lo)]  hy  using  Cayley’s  method;"® 

^  l±hQ^ 

1  —  ^iQAx 

If  the  Q  operator  is  now  discretized  by  a  finite-difference  approximation  in  i ,  a  complex 
tridiagonal  system  of  equations  results.  This  is  simply  the  standard  Crank-Nicholson 
method  for  a  parabolic  partied  differential  equation.^' 

There  are  two  fundamental  problems  with  this  approach  however.  First,  the  trmica- 
tion  error  of  the  Crank-Nicholson  scheme  is  only  second  order  in  Ac  and  Ax,  and  thus 
reqxiires  small  mesh  intervals  to  achieve  iiigh  accuracy.  This  in  turn  leads  to  very  large 
matrix  systems  and  many  range  steps.  The  second  problem  with  the  finite-difference  ap¬ 
proach  is  a  lack  of  rigorous  error  bounds  on  the  solution  that  can  be  monitored  during  the 
comptitations  to  asstne  a  fixed,  preset  global  error. 

To  deal  with  both  these  issues,  an  alternative  solution  method  is  tised  that  is  based  upon 
spectral  operator  techniques.  This  spectral  technique  is  known  as  the  split -step  Fourier 
PE  algorithm  (SSFPE)  and  w'as  developed  by  Hardin  and  Tappert.^  The  remainder  of 
this  section  deals  with  development  of  the  SSFPE  algorithm,  associated  error  bounds,  and 
numerical  implementation  using  the  fast  Fourier  transform  algorithm. 

3.1  Magnus  Expansion 

Given  the  parabolic  equation  Eq.  (78),  there  still  exists  a  problem  of  developing  a 
numerical  solution.  The  difficulty  arises  from  the  nonlocal  nature  of  the  Q  operator.  To 
solve  the  GPE,  techniques  from  time-dependent  quantum  scattering  theory  are  used  to 
express  the  solution  in  terms  of  an  evolution  operator.*^®  Let  u.s  define  the  wave  evolution 
operator  C/(x,xo)  that  determines  the  PE  wave  field  ri(x)  at  the  point  x  =  (x,  z)  in  terms 
of  the  wave  field  ri(xo)  at  the  point  xo  ==  (j’Ot'^)  by  the  equation 

^(x)  =  r(x,xo)0{xo). 

It  follows  by  substitution  into  Eq.  (78)  that  U  is  an  operator  satisfying  the  partial  differ¬ 
ential  equation 

^(/(x.xn)  ....  .  _ 

Clearly,  U  must  also  satisfy  the  initial  condition 

fTxO’Xo)  =  1.  (85) 
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and  possess  the  group  property 


y(x,X(i)  =  P(x,xi)f/(xi,xo). 


(86) 


Forthermore,  if  there  are  no  dissipative  processes  present,  it  follows  from  the  complex 
Poynting  theorem  that  the  ^-integrated  field  energy  density  must  be  constant  with  range 

J  j0(a;,?)pd2  =  J  \xl){xo,z)\^dz, 

and  this  leads  to  the  important  result  that  the  evolution  operator  U  must  be  unitarj^: 
I7^(x, Xo)  =  C/~^(x,xo),  where  is  the  (Hermitian)  adjoint  of  U.  The  importance  of 
having  a  unitary  operator  for  numerical  work  should  not  be  underestimated,  since  this 
form  for  U  leads  to  numerical  methods  that  can  be  shown  to  be  absolutely  convergent  for 
all  range  step  sizes. 

Integrating  Eq.  (84)  and  using  Eq.  (85)  yields  an  integral  equation  for  the  evolution 
operator 

?7(x,xo)  =  l  +  t  /  Q{t,z)U(t,XQ,z)dt, 

JXQ 

Solving  the  above  integral  equation  iteratively  yields  Dyson’s  expansion  for  U: 


fz  fx  rx-i 

U(x,XQ,z)=il  +  i  I  Q(t,z)dt-I  dxi  I  dx2Qixi^z)Qyx2,z) 

Jxq  J-:o  dxo 

fX  fXl  fX2 

-i  I  dxi  I  dx2  I  dx3Q{xi,z)Q{x2,z)Q{xz,z)  + 
JxQ  Jxq  Jxo 


(87) 


The  x-ordering  of  the  integrands  in  Eq.  (87)  is  very  important  because,  in  general,  the 
Q  operators  at  different  x-positions  do  not  commute.  Thus  the  ordering  of  the  operators 
is  determined  by  the  range  to  which  they  refer,  with  operators  referring  to  earlier  ranges 
always  appearing  to  the  right.  Though  the  Dyson  expansion  is  a  formally  exact  solution 
for  U,  it  is  not  very  useful  for  calculations,  since  if  the  series  in  Eq.  (87)  is  truncated  after 
a  finite  number  of  terms,  the  result  is  no  longer  a  unitary  operator. 

A  formal  solution  to  Eq.  (84)  that  does  preserve  the  unitarj  nature  of  U  can  be  found  by 
using  a  technique  developed  by  Magnus.^^  The  Magnus  result  is  based  upon  en  exponential 
operator  expansion  which  effectively  includes  all  terms  in  the  Dyson  expansion.  The  details 
are  found  in  Ryan,®  and  yield 


U(x,XQ,z)  -  expfi  f  Q{xi,z)dxi  f  dxi  f  dx2[Q{xi,z),Q(x2,z)] 

\  Jxo  ^  Jxq  Jxo 

^  fX  rxi  rx2  f 

-hi  dx^  dx2  I  dx3  <  [[Q(xi,2),Q(x2,2)],Q(x3,z)]  (88) 
O  Jxn  Jxo  Jxo  V 

+  [[<3(X3,2),Q(X2,^)],(?(X1,2)]|  4-  •••  Y 
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While  the  Magnus  expansion  provides  an  exact  formal  solution  to  the  GPE  wave  field, 
three  approximations  must  be  made  before  it  is  useful  for  numerical  computations.  First, 
the  infinite  series  in  the  Magnus  exponent  must  be  truncated.  Keeping  the  first  term  in 
the  operator  exponential  in  Eq.  (88),  yields  the  first-order  Magnus  approximation  Ui : 


U{x,xo,z)‘:^Uiix,XQ,z), 

Ui{x,XQ,z)  =  =  ^iix-xo)H^ 


where  H  is  the  range-averaged  PE  “Hamiltonian” 

j  rxQ+/S. 

H  = I  Q{t,z)dt,  A  —  x  —  xq. 
^  Jxq 


(89) 

(90) 


Second,  the  pseudodifferential  GPE  operator,  Q,  is  factored  into  the  sum  of  two  ordinary 
operators,  one  of  which  is  independent  of  range; 


Q(x,z)  Rs  .4(2)  -P  B{x,z). 

(91) 

This  factorization  of  Q  is  not  unique;  in  general,  the  A  and  B  operators  do  not  commute 
with  each  other.  For  example,  if  the  SPE  operator,  Eq.  (82),  is  used  to  approximate  Q{x), 
then 

Ai-)-  ^ 

2hdz^^ 

Bix,z)^{ko-k)-^V{x,z), 

(92) 

(93) 

while  if  the  WAPE  operator,  Eq.  (83),  is  tised 

A{z)  -  \/kl  +  ^~ko 

(94) 

B{x,  z)  =  K{x,  z)  —  k. 

(95) 

Finally,  the  Hamiltonian  H  is  evaluated  by  using  the  midpoint  rule  as 

PEQ-fA 


1  fXo+i\ 

Hix,z)==--  Q{t,z)(H, 

^  Jxo 


Q(xo  +  A/2,  z)  -p  ,  t  €  (xQ,  xo  +  A), 

A{z)-^B{x,z)  x  =  xo-i-A/2. 


(96) 


If  this  is  done,  a  formal  solution  of  Eq.  (78)  can  be  expressed  in  exponential  operatoi 
form  as 

ip{x  -h  A,^)  ~  Ui{x  -P  A,x,z)‘il){x,z), 

_  g-Pi AH(x,  z)  2), 

r 

=  exp  |■P^[AA(^:)  -p 
Rs  exp  {-l-iA{A(z) -P  B(x.  z)]}  V’(x,r).  (97) 


r 


Bit,z)dt]  >  t/>(x,;), 
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3.2  Split-Step  PE  Algorithm 


The  first-order  Magnus  expansion,  Eq.  (97),  -s  not  suitable  for  numerical  work  since  the 
A  and  B  operators  appearing  in  the  exponent  do  not  in  general  commute.  To  resolve  this, 
the  Trotter  product  formula^®  is  used  to  symmetrically  factor  the  Magnus  expansion  into 
the  product  of  simpler  operators,  C/^  and  Vq 

Ua  =  i^.4(A)  = 

Ub  =  Ub{A)  = 

to  yield  the  split-step  PE  algorithm: 
tl}(x  +  A,2)  ~ 

_  g-htAA(.-)/2  g-»-iAS(a+A/2,2)  g-f-iA.4{3)/2 

The  operator  17^  is  known  as  the  free-space  propagator  and  is  the  exact  formal  solution  to 
the  parabolic  equation  in  the  absence  of  refraction.  Physically,  the  split-step  PE  algorithm 
amounts  to 

(1)  a  half-step  of  free-space  propagation, 

(2)  a  phase  correction,  l/s(A),  to  account  for  refractive  effects;  and 

(3)  finally,  another  half-step  of  free-spEice  propagation,  Ua{y)- 

In  actual  use,  the  split-step  algorithm  is  used  to  advance  the  PE  field  multiple  range 
steps: 

tpn  =  e*An^/2  g»A„S„  An  A/2  ^ 

where 

An  —  Xfi  —  Xjn—Xi 

t/’n  ~ 

Bn  -  B[(a:n-1  +  3:n)/2,  z] . 

By  using  the  group  property  of  the  evolution  operator,  Eq.  (86),  the  split-step  algorithm 
may  be  iterated  to  give 

■A/2gt'A5n  g*  An  4/2^>  A,,  — 1  i4/2^»  A5n  —  1  An  — 1  A/2  2 

_  giA„  A/2^«ABng*{'^n-f  A„_i).4/2^iAB„_igjAn-i>4/2 

_  ^tA„4/2gtAnSng‘(An+An-i)A/2  , ,  _giA2B2gi(A2+Ai)4/2giAifiigiAi^/2 

3.3  Split'Step  PE  Truncation  Error 

At  each  step  of  the  PE  calculation,  there  will  be  a  certain  amount  of  error  induced  in 
the  solution  due  to  round-off  error  and  the  intrinsic  error  of  the  split-step  PE  algorithm. 
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The  round-ofF  error  is  caused  by  performing  calculations  using  floating-point  arithmetic  on 
finite- word" length  computers.  This  type  of  error  is  highly  dependent  upon  the  hardware 
characteristics  and  actual  numerical  algorithms  used,  and  will  not  be  discussed  further. 
The  intrinsic  error  is  termed  truncation  error  and  arises  from  the  various  approximations 
in  going  from  the  exact  Magnus  solution  in  Eq.  (88)  to  the  split-step  algorithm  in  Eq.  (98). 
The  split-step  PE  is  thus  predicated  on  four  approximations: 

(1)  approximating  the  exact  square  root  GPE  operator,  Eq.  (79),  by  a  finite  sum  of 
ordinary  operators; 

(2)  truncating  the  fonnally  exact  Magnus  solution  for  the  evolution  operator  and  keep¬ 
ing  only  the  first  term  in  the  exponent,  Eq.  (89); 

(3)  evaluating  the  PE  Hamiltonian  by  the  midpoint  rule,  Eq.  (96);  and 

(4)  approximating  the  exponentiM  Magnus  operator  by  a  symmetrized  factorization  of 
individual  operators,  Eq.  (98). 

Each  of  these  steps  introduces  an  intrinsic  truncation  error  in  the  split -step  PE  solution — 
the  abihty  to  quantify  these  errors  is  crucial  to  the  successful  implementation  of  a  numerical 
PE  code. 

Although  the  Magnus  expansion  provides  an  exact  formal  solution  to  the  PE  wave  field, 
three  approximations  must  be  made  to  it  before  numerical  computations  are  feasible.  First, 
the  infinite  series  in  the  Magnus  exponent  is  trancated  after  the  first  term 

t/(xo  +  ^  -  T  I  Q{t,z)dt,  (99) 

^  Jxo 

where  H  is  the  range- averaged  PE  Hamiltonian. 

The  truncation  errors  incurred  by  using  Eq.  (99)  are  a  function  of  the  range  dependence 
of  the  Q  operator.  If  Q,  or  equivalently  AT,  is  independent  of  x,  then  Eq.  (99)  is  exact.  For 
range-dependent  K,  the  truncation  error  in  U  may  be  quantified  by  examining  the  next 
term  in  the  Magnus  expansion  and  gives 

{XO  +  A  ^  "I 

f A/f  -  i  J  dx\  j  d2-2[(5(.'ri,2),Q(x2^^)]  /  . 

=  exp  ^  •  (100) 

Thus,  when  Q  is  range-dependent,  truncation  errors  proportional  to  O(A^)  are  incurred 
no  matter  what  form  of  the  PE  operator  is  used.  Furthermore,  how  the  range-averaged 
PE  Hamiltoniein  H  is  evaluated  critically  affects  the  error  budget.  For  example,  if  Q  is 
expanded  in  a  Taylor  series  about  the  midpoint,  x  =  {x+xq)/2,  and  the  integral  in  Eq.  (98) 
is  evaluated  by  the  midpoint  rule,  then 

U{xo  +  A,a:o,2)  =  exp(i/\Q{x,z)  +  ~  {iQ"(i,z)  -  4[Q'{x,  z),Q{x,  z)]  }  -f-  0(A‘^)^, 
where  Q'{x,z)  —  dQ{x,z)ldx  and  Q"{x,z)  =  d^Q{x,z)/dx‘^ . 
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The  local  truncation  error  Ti  incurred  in  the  PE  field  solution  by  using  Eq.  (99),  with 
Q  evaluated  at  the  midpoint  x,  is  then 

Ti  =  +  A,xq,z)  ~  Uxo^z), 

A  3  , 

-  z)~4:  [Q'(x,  z)Q{x, z)  -  Q(x,  z)Q'{x,  2)]  j?/)(xo,  z) , 
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{' 


dx^ 


— -  ^(^)— g— j  I  ^(xo,  z) . 


(101) 


A  second  type  of  truncation  error  occurs  if  the  symmetrized  operator  splitting,  defined 
by  Eq.  (98),  is  employed  to  compute  the  propagator.  This  splitting  error  can  be  evaluated 
by  using  the  Baker-Campbell-HausdorfF^^~^^  (BCH)  expansion  of  two  noncommuting 
operators;^^ 

gjtgB  _  4.  5  q.  |[.4,S]  +  ]^[A,  [A,B]]  +  ■^[[A,  B],B]  +  •••). 

Applying  the  BCH  expansion  to  Eq.  (98)  gives 

=expi|AA  +  AB-  +  OCA-*)!  , 


and  the  local  tnmcation  error  T2  caused  by  using  the  Trotter  product  formula,  with  B 
evaluated  at  the  midpoint  is  thus 


72  = 


;A[24+B(*)]  _  p»AA/2  giAB  ^iAA/2 


1a 

iA^ 

'W 


{2[[A,B(x)],i?(i)]  -  [[B{x),A],A]]^k{xo.z), 

[2(AB2  _  2BAB  +  B^A)  -  (FA^  -  2ABA  +  A^B)]  0(xo,  z) . 


(102) 


If  the  SPE  operator  is  used  to  approximate  Q,  then  the  third-order  local  truncation 
error  terms  are  given  explicitly  by 


A^  d^V  _  d^V 


-  48  (**°5x2  '^^dz'^dx 


^  0(XO,2), 


(103) 


and 


T  - 

^  mo 


l^V 

4  d. 


d^v  ,2(^yy 

z^’^'^^dz’^dx  ^^\dz) 

ei^v  df^v 


+ 


i2ko 


dzdx  dz^ 


V’(xo,z) 

dil^{xo,z)  d’^V  d^iixo^z) 


dz 


+ 


522  5^2 


}■ 


(104) 
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where  V  —  V(x,z).  Equatioiis  (103)  and  (104)  show  that  the  local  truncation  error  for 
the  symmetrized  operator  splitting  is  linear  in  frequency  (via  the  kfj  terms)  and  cubic  in 
the  range  step-size,  A. 

To  simplify  the  an2ilysis,  let  us  assume  that  the  potential  energy  function  V  is  inde¬ 
pendent  of  the  range  x.  (In  most  applications,  the  horizontal  ^-derivatives  of  V  are  much 
smaller  than  corresponding  vertical  z-derivatives.)  In  this  case,  the  third-order  local  ab¬ 
solute  error  in  the  PE  field,  caused  by  advancing  rvith  the  range  step  A  is  given 

by 


S^^^'ip{xo  -h  A,2)  = 


iA^  (\ld^v  ,2fdV\^ 
48A:o\[4a24  “VOrJ 


4ixo,z) 


d^Vd'4ixo,z)  ,  d^V  d^7j>{xo,z)\ 


Instead  of  working  with  the  pointwise  error  estimate,  defined  above,  it  proves 

more  useful  to  consider  the  totcil  relative  error  per  range  step  which  is  normalized 

by  the  vertical  energy  density  of  the  PE  field: 


p(3)V>il 


__  f  i/)* tp  dz 
Jli'l'^dz 

The  term  may  be  evaluated  by  integrating  by  parts  and  takes  the  form 


(106) 


<^vdtp 


+ 


dvy 


1  d'^V 

a,  i  - 

^  dz^ 


dzj  dz^ 

aiV’(a:o,0)l^ 


dz 


dz^  dz  dz^  dz^ 
2 


dz , 


1  d'^V 

4  5^2  dz^ 


dz 


r=0 


dz 


}• 


(107) 


To  ensure  that  the  total  error  in  the  computed  PE  solution  is  bounded,  the  PE  step-size 
A  is  chosen  so  that  the  third-order  error  term  is  small  relative  to  the  first-order  term 

£;(3) 

^(iT  -  ^  ’ 

where  t,.^i  is  a  specified  error  tolerance,  and 

M)  _  MllVll 
IK'U  ’ 


d'^ 


4'  dz . 


(108) 
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This  implies  that  the  PE  range  step-size  should  be  chosen  so  that 


<  24erei 


and  determines  the  PE  range  step,  A,  by 

I  V 
A2 


24e„,  Jr  -  iiV 


^  dz 


J 


kl  (^\ 

“UJ 


d4)  ^ 

dz 

^dz^  dz^ 

^4  3z2  dz 

z-Q 

(109) 


In  numerical  implementation,  the  vertical  derivatives  appearing  in  Elq.  ( 109)  are  evaluated 
by  numerical  finite-difference  approximations.  As  the  PE  code  advances  the  field,  the  local 
error  budget  is  monitored  and  the  range  step-size  is  dynamically  adjusted  to  keep  the  local 
error  below  a  preset  threshold.  The  question  of  a  more  detailed  error  analysis  will  not 
be  covered  in  this  report.  The  interested  reader  is  referred  to  Ryan*  for  analysis  of  the 
split-step  PE  truncation  errors. 


3.4  Split-Step  Fourier  PE  Algorithm 


The  PE  solution,  A,  2),  at  range  a; -FA  is  obtained  from  the  known  field,  0{x,  z),  at 
range  x  by  means  of  the  split-step  PE  algorithm,  Eq.  (98).  The  presence  of  the  differential 
operator,  A(z),  in  the  exponent  can  be  dealt  with  by  transforming  to  a  basis  in  which 
A{z)  is  diagonal.  One  such  basis  is  the  Fourier  basis,  with  the  z-space  Fourier  transform 
T  being  defined  as 


/oo 

dz, 

•cx> 


and  the  inverse  transfomi  T  ^  being  defined  by 

1 


The  conjugate  transform  variable  p  can  be  associated  with  the  vertical  wave  number,  via 
p  =  ibosin^,  where  B  is  the  propagation  angle  with  respect  to  the  horizontal. 

As  is  well  known,  the  vertical  derivative  operator  in  z-space  is  related  to  the  conjugate 
p-space  operator  by 


(-?p)", 


in  which  case  the  SPE  operator  Eq.  (S2)  is  implemented  via 


(110) 
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Similarly,  the  WAPE  operator  Eq.  (94)  is  implemented  via 


-if  z)] 


=  > 


(111) 


Now  any  function  may  be  expressed  in  terms  of  even  and  odd  components,  so  the  PE 
field  can  also  be  written  as 

^(l,  z)  -  tj)e{x,  z)  +  ^0(3:,  2), 


where 


^c(ar,  z)  ^  z)  +  tk{x,  -2)], 

0o(a?,  z)  =  ll^ix,  z)  -  xl}(x,  -z)]. 

Being  a  linear  operator,  the  Fourier  transform  may  also  be  expressed  in  terms  of  even  and 
odd  components  as 

4'(x,p)  =  T[ii>{x,z)]  =  +  ’J'o(a;,p), 

where  ami  are  represented  in  terms  of  Fotirier  sine  and  cosine  transforms  as 

roo 

'^e{x,p)  =  fc{'tl’e{x,z)]  =  2  /  t{?eix,z)  COS pzdz, 

Ja 

and 

roo 

vpo(x,p)  =  ^a[^o(x,2)]  =  ~i2  I  xl)o{x,z)  sin pzdz. 

Jo 

The  corresponding  inverse  sine  and  cosine  transforms  are  defined  by 

1  /"°® 

ijuix,z)  =  jF"^[^e(a^,p)]  =  -  /  ^e(a;,p)cosp2<ip, 

TT  Jo 

and 

i 

i^o{x,z)  =  i^7^[^o(a:,p)]  =  -  /  j{x,p)  sin pz  dp. 


In  practice,  the  infinite  Fourier  transforms  are  replaced  by  finite  discrete  Fourier  trans¬ 
forms  (DFTs)  over  the  domain  (0,  Zmax),  which  in  turn  are  evaluated  numerically  using  the 
fast  Fourier  transform  algorithm.  The  actual  FFT  algorithms  implement  fast  real- valued 
sine  and  cosine  transforms  to  reduce  core  storage  and  improve  computational  speed. 

When  using  DFTs,  it  is  important  to  satisfy  the  Nyquist  sampling  theorem  to  avoid 
transform  aliasing  problems.^^  If  Pmax  is  the  maximum  vertical  wave  number  and  N  is 
the  discrete  cosine/sine  transform  size,  then  the  Nyquist  condition  is  just 

Zmax  Pmax  —  •  [112) 


30 


3.5  Reflectioniess  Absorber  Boundary 


The  appropriate  boundary  condition  to  be  satisfied  as  ^  -4  oo  by  the  PE  field  ip(x,  z)  is 
the  Sommerfeld  outgoing  wave  radiation  condition,  Eq,  (35).  Since  the  split-step  Fourier 
algorithm  employs  finite  Fourier  transforms,  the  implementation  of  a  radiation-type  bound¬ 
ary  condition  is  quite  complicated.  This  follows  from  the  fact  that  truncation  of  the  infinite 
^-domain  down  to  a  finite  interval  in  the  Fourier  transform  leads  to  the  introductioi-  of 
spurious  discrete  standing  wave  solutions  in  the  vertical.  In  effect,  the  terminal  impedance 
at  the  end  of  the  transform  grid  is  not  properly  “matched”  to  the  radiation  boundary 
condition. 

To  circumvent  this  problem  and  attenuate  the  spurious  standing  wave  solutions  intro¬ 
duced  by  the  finite  Fourier  transforms,  a  complex  absorber  potential  Vahsi^)  is  added  to 
the  split-step  B  operator: 


B{x,z)  =»  B(x,z)  +  —Vabsi^)- 

The  particular  form  of  the  complex  absorber  or  “sponge”  is  found  by  recourse  to  a  nuclear 
optical  model^®  analog  based  upon  the  theory  of  reflectionless  potentials,^®  The  specific 
form  chosen  is 

^"a6a(^}  =  V^0Sech2[a(2 

^max  )]  (113) 

where  the  parameters  {Vh,a}  are  determined  parametrically  by  minimizing  transmission 
and  reflection  coefficients  from  the  sponge  region.'*® 

3.6  PE  Starting  Fields 

The  split-step  PE  method  must  be  initialized  with  a  starting  field  distribution 
at  some  distance  iq  from  the  source,  since  the  parabolic  wave  equation  is  not  valid  in  the 
region  of  the  source.  Two  options  are  available  to  compute  the  starting  fields.  First,  the 
initial  «-space  field  may  be  obtained  by  analytic  methods,  assuming  free-space  conditions 
and  treating  the  source  as  a  point  radiator.  This  option,  though,  is  not  useful  for  highly 
directional  antennas. 

The  second  approach  is  to  use  the  duality  of  the  antenna  aperture  field  distribution  and 
antenna  radiation  pattern  function.**  In  free  space,  the  antenna  radiation  pattern  function 
f{p)  and  the  antenna  aperture  field  distribution  A{z)  are  a  Fourier  transform  pair: 

/OO 

dz ,  (114) 

-OO 


where 


p=  ko  sin  6, 

pO  =  fco  sin^o- 
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(115) 

(116) 


In  Eq.  (115),  Q  is  the  elevation  angle  measured  with  respect  to  the  horizontal  (6  >  0  is 
up),  and  in  Eq.  (116),  is  the  antenna  mainlobe  vertical  pointing  direction.  The  antenna 
pattern  firnction  is  used  to  construct  even  and  odd  symmetry  p-space  fields  in  the  form 

S'eiO.p)  =  +  /(-?)«+■'"“ ,  (117) 

and 

’5o(0,p)  -  /(p)e-*P‘-'>  -  /(-p)e+*>° ,  (118) 

where  the  Fourier  shift  theorem  has  been  used  to  properly  account  for  a  nonzero  source 
altitude,  zq.  Given  the  above  even/odd  symmetry  p-space  fields,  the  corresponding  2-space 
PE  field  is  obtained  by  taking  the  inverse  cosine  or  sine  transform  of  Eq.  ( 1 1 7)  or  Eq.  (118), 
respectively: 

^>e(0,2)  =  :F7^[$e(0,p)], 


and 


t/ao,-)=57M^o(o,p)]. 

3,7  Antenna  Patterns 

The  PE  method  is  capable  of  modeling  the  radiation  emitted  by  directional  antennas, 
provided  that  the  complex  antenna  radiation  pattern,  Eq.  (114),  is  known.  Often  this 
information  is  not  av,  ilable  for  specific  radar  systems,  so  simplified  generic  antenna  pat¬ 
terns  are  used.  These  generic  patterns  display  some  of  the  features  of  real  antennas  while 
retaining  fairly  simple  analytical  forms.  Each  of  the  analytic  antenna  patterns  is  specified 
in  terms  of  a  normalized  p-space  steering  parameter,  t,  defined  by 

^  _  sin  9  —  sin  _  P  ~  PO 
sin(5^ixi)  Pi 

where  ^ 

Pi  =  ^’osjn(o^6w')i 
2  Z 

and  is  the  half-power  (i.e.,  3-dB  down)  beamwidth.  The  follow'ing  analytic  antenna 
radiation  patterns  are  useful: 

3.7.1  sin(x)/x 

The  imiformly  illuminated  aperture  corresponds  to  a  radiation  pattern  having  the  func¬ 
tional  form 

f{t)  =  0  =  1.3916.  (119) 

at 

The  scale  parameter  a  is  determined  by  solving  the  nonlinear  equation 

sinr  1 

T”  “ 

This  pattern  has  the  narrowest  mainlobe  width,  at  the  expense  of  high  sidelobe  levels. 
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3.7.2  Gaussian 


The  Gaussian  antenna  radiation  pattern  has  the  functional  form 

f(t)  ~  e~^  ^  ,  where  (120) 

This  pattern  is  the  optimal  compromise  between  sidelobe  level  and  mainlobe  width 

3.7.3  Compound 


The  compotind  radiation  pattern  is  a  one-parameter  pattern  that  is  a  blend  of  the 
uniformly  illuminated  aperture  and  the  cosine-squared  aperture  distribution,  having  the 
functional  form 


m== 


2  sm(at) 
1  +  c  at 


1-c  1 

2  1  — (ai/7r)2 


0  <  c  <  1. 


The  parameter  a  controls  the  mainlobe  width  and  is  determined  by  solving  the  nonlinear 
equation 


sin  a 
a 


2c  + 


l-(a/?r)2j 


1  +  c 

V2  • 


The  uniform  aperture  corresponds  to  c  =  1,  while  the  cosine-squared  aperture  corresponds 
to  c  =  0. 


3.7.4  Hansen 

Another  single  parameter  pattern  is  obtained  from  the  circular  aperture  distributions 
analyzed  by  Heinsen.^^  This  pattern  has  the  fimctioual  form 


/w  = 


H 

ii{H) 


h(V») 
^  ’ 

ii(vlid) 

y/F\ 


for  a*  >  0 
for  a-  <  0, 


(121) 


where 

X  =  —  (ot)^, 

and  ji  and  ii  are  the  first-order  spherical  and  modified  Bessel  functions  defined  by 

sin  z  cos  2 


.  ,  ,  sinhr  cosh* 

*l(2)  = - ^  + - • 

z 

The  quantity  a,  which  determines  the  3-dB  point  in  the  pattern,  is  found  by  solving  the 
transcendental  equation 


\/2H 


for  C  >  j 
for  C  <  5  . 
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The  first  sidelobe  level,  SLL,  in  the  radiation  pattern  caoa  be  shown  to  have  a  value  of 


30.84  +  20  log 


'3ii(iiy 

H 


dB, 


down  from  the  peak,  and  is  located  at 


The  parameter  H  allows  a  trade-off  between  low  sidelobe  level  and  mainiobe  beamwidth. 
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§4.0  CONCLUSION 


This  report  d^ribes  the  basic  physics  2uid  analytical  techniques  used  to  model  mi¬ 
crowave  propagation  in  range-dependent  environments.  These  techniques  have  been  im¬ 
plemented  in  a  range-dependent  propagation  computer  model  —  the  VTRPE  code.  The 
VTRPE  computer  model  is  used  to  predict  and  analyze  the  performance  of  radar  and  com¬ 
munication  systems  operating  in  spatially  varying  environments.  The  methods  described 
herein  account  for  spatial  variations  in  the  atmospheric  refract  ivity  as  well  as  for  variations 
in  surface  dielectric  properties  and  topography. 
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